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ABSTRACT 

We demonstrate how the Syer & Tremaine made-to-measure method of stellar dy- 
namical modelling can be adapted to model a rotating galactic bar. We validate our 
made-to-measure changes using observations constructed from the existing Shen et al. 
(2010) N-body model of the Milky Way bar, together with kinematic observations of 
the Milky Way bulge and bar taken by the Bulge Radial Velocity Assay (BRAVA). 
Our results for a combined determination of the bar angle and bar pattern speed 
(~ 30° and ps 40 km/s/kpc) are consistent with those determined by the N-body 
model. Whilst the made-to-measure techniques we have developed are applied using 
a particular N-body model and observational data set, they are in fact general and 
could be applied to other Milky Way modelling scenarios utilising different N-body 
models and data sets. Additionally, we use the exercise as a vehicle for illustrating 
how N-body and made-to-measure methods might be combined into a more effective 
method. 

Key words: Galaxy: kinematics and dynamics - Galaxy: structure - methods: N- 
body simulations - methods: numerical 
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1 INTRODUCTION 

The made-to-m e asure method (M2M) proposed by 
ISver fc Tremainel |l996) for modelling stellar dynamical 
systems has been the subject of growing interest recently 
(for example, Jpurdcuil & Emscllcm 2007. de Lorenzi et al 



2007 



2010 



de Lorenzi et al l |200S|, iDehnenl 12009] lLong fc Mac 



Das et all l201ll . iMorganti fc Gerhard! 



20121). This 



paper is the third in a s er ies to gether with lLong fc Mad 
(|2010l ) and lLong fc Mad" (|2012t ). Whereas the previous 
two papers, in common with other authors, dealt with 
theoretical galaxy models and external gala xies, this paper 
applie s the M2M method to the Milky Way. iBissantz et al.l 
|2004 ) previously used the method to model the equatorial 
plane surface brightness of the Galaxy. Note however 
that the effectiveness of their model is subject to debate 
|Rattenburv et ail [2007). We believe that the current 
investigation is the first in which Galactic M2M modelling 
has taken place using a rotating frame, and with luminosity 
and kinematic observables based on fields specified in terms 
of Galactic longitude and latitude. 



IShen et all (^Old - ) investigated whether the Milky Way 
bulge mig ht be formed from disk and bar related mecha- 
nisms. |Shen et al.l I20l61 provide a much fuller analysis of 
the issues.) Starting from a simple N-body disk simulation of 
some 10 6 self-gravitating particles, a bar develops naturally 
(from the conventional bar instability), buckles and thick- 
ens, and a pseudo-bulge appears. With appropriate scal- 
ing, they are able to match kinematic data take n from the 
BRAVA (Bulge Radial Velocity Assay) surv ey l|Rich et al.l 
l2007l : lHoward et al.|[20oi : lKunder et al.ll2012r i. and thus sug- 
gest that the formation of the Milky Way bulge is mainly 
shaped by the internal dynamical buckling instability. Their 
model weakly constrains the Galactic bar to an angle of 
m 20° to the Sun-Galactic centre line, with a pattern speed 
(J. Shen private communication ) of ~ 40 km/s/kp c. In our 
M2M modelling, we utilise the IShen et all i|2010h particle 
data and the BRAVA field observations, and investigate 
whether the values for the ba r angle and patte rn speed we 
determine are consistent with lShen et al.1 (|20ld h Note that 
our intention is not to refine or improve upon the N-body 
results. 



* E-mail: rjl2007@gmail.com 



The M2M method has been demonstrated to be a 
flexible tool for modelling with a variety of observables 
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with measurements distributed in various different spa- 
tial configurations. The method can handle luminous and 
kinematic observables and there is no requirement that 
the observable measurements should share a common spa- 
tial positioning. However in order to function it does re- 
quire a system of particles with a set of initial conditions 
and a gravitational potential to generate their orbits. N- 
body models have been used extensively to model galactic 
component r y, for example disks, spirals, b ars, black holes 
llHernauistl ll993l, ISigurdsson et all Il995l. ISeUwoodl |2012| . 



ISellwood fc Wilkinson! Il993l . iFuxl Il999l . IShen et"aiTl2010l ) 
The models are effective in modelling the general proper- 
ties of generic galaxies but have limited tailoring capabilities 
for modelling the kinematics of real galaxies. Thus, there 
are potentially advantages to be gained by combining the 
two modelling methods with the N-body method providing 
an initial model including particles and potential, and the 
M2M method refining the model by weighting the particles 
to match the measured observables. One might imagine a 
sequence of interleaved N-body and M2M models with feed- 
back from one model type being used to improve modelling 
with the other. By checking consistency or otherwise of re- 
sults between the two approaches, the current paper is an 
initial step in this direction. Note that we will not be utilising 
interleaving as described earlier but using an existing pub- 
lished N-body model and seeking to use the M2M method 
to corroborate the N-body model and results. If corrobora- 
tion can not be achieved for this simple case then we must 
understand (and address) the underlying reasons before de- 
veloping a more tightly coupled approach between the two 
methods. 

To summarise, our objectives in performing this inves- 
tigation are three-fold: 

(i) to update our M2M implementation, and the theory 
underlying the M2M method if necessary, to support mod- 
elling of a rotating galactic bar and to handle Milky Way 
observations, 

(ii) to validate our M2M changes using the [gh en et alj 
i|2010t ) Milky Way N-body bar model and the BRAVA obser- 
vations with the expectation that we should achieve results 
consistent with the N-body model, and 

(iii) to aid our understanding of the synergy that might 
be developed between coupled N-body and M2M modelling. 

The structure of the paper is as follows. In section [5] we 
describe the M2M method, and in section [3] we consider 
the changes required to construct and run our M2M bar 
models. In sections 2] and HI] we discuss our results and draw 
conclusions within the context of the objectives. In section[5] 
we discuss further combining N-body and M2M modelling. 



2 THE M2M METHOD 
2.1 Outline 

In brief, the M2M method is concerned with modelling stel- 
lar systems and individual galaxies as a system of test par- 
ticles orbiting in a gravitational potential. Weights are as- 
sociated with the particles and are adjusted as the particles 
are orbited so that, by using these weights, observational 
measurements of a real galaxy are reproduced. As conceived 



an d demonstrated bv l Sver fc Tremaind (|l996l ) and extended 
bv lLong fc Mad (|2010h (quantifying weight convergence) , we 
expect that the weights themselves will have converged in- 
dividually to some constant value. It is natural to relate the 
particle weights to the luminosity of a galaxy and then to 
consider how the galaxy's surface brightness and luminosity 
weighted kinematics could be generated using the particle 

system. 

In the next section, based heavily on lLong fc Mad 

we set out the theory underlying the M2M method. 



2.2 Theory 

For a system of TV particles, orbiting in a gravitational po- 
tential, with weights Wi, the key equation which leads to the 
weight evolution equation is 



(i) 



where % , S and d are all functions of the particle weights 
w = (ttfi, • • -,tojv); t is time; and \i and e are positive 
parameters. Q is the number of additional constraints d. 
The equations governing weight evolution over time come 
from maximising F(w) with respect to the particle weights 
(dF/dwi = 0, Vi) and rearranging terms to give equations 
of the form 

^w. = -£WiG(w). (2) 

The overall rate of weight evolution is controlled by e. The 
precise form of the function G(w) depends on the constraints 
d and is illustrated later (equation I10|) . The process be- 
ing applied to \ 2 i s one °f regularised, parameterised con- 
strained extremisation. 

The \ 2 term in F arises from assuming that the proba- 
bility of the model reproducing a single observation can be 
represented by a Gaussian distribution and then construct- 
ing a log likelihood function covering all observations. For 
K different observables, we take \ 2 in the form 



X 



(3) 



where X k are small, positive parameters whose role is ex- 
plained in section [3~6l Note that for statistical purposes, for 
example obtaining confidence limits or regions, equation [3] 
means that \ 2 is a linear combination of \ 2 terms. As such, 
it must not in general be associated with a \ 2 distribution 
but with a sum of gamma distributions. 



2 

Xk 



3 



and 



y k ,j{w) - Y kJ 



(4) 



(•») 



where Y k ,] is the measured value of observable k at position 
j with error a k ,j , and yt,j (w) is the model equivalent of Yu ,j ■ 



(w) = WjKk,j(ri, Vj)8(i G k,j) 



(6) 
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where Kkj(ri,Vi) is the kernel for observable k evaluated 
at position j for a particle with position n and velocity Vi. 
8(i € k, j) is a selection function and signifies that only par- 
ticles which contribute to observable k at position j should 
be included in the calculation of yuj- We have listed the 
kernels required for this paper in section ET21 
The entropy function S in F is 



JV 

S(w) = -\^ u , i l n (B.) 



(7) 



where rrii is taken as the initial value of a particle weight (in 
practice, we take rrii = 1/N). S is used for regularisation / 
smoothing purposes with the amount of regularisation being 
controlled by the parameter /x. The derivative term dS/dt 
indicates that over time we require the particle weights, and 
thus 5, to be constant. As demonstrated in lLong fc Mad 
(2010), the term behaves as the constraint dS/dt = 0. The 
behaviour of the regularisation terms in the weight evolu- 
tion equations i s to move the particle weigh ts towards rrii /e 
rather than m,. lMorganti fc Gerhard] |2012l ) solve this prob- 
lem by changing S to 



S(w) = 



N 



L rrii 



- 1 



(8) 



Modifying 5* in this way causes no change to equation [T] 
Other regularisation schemes are possible, for example sec- 
ond derivative smoothing in integral space. The challenge 
then is to produce a cost-effective implementation in a mes- 
sage passing computational environment. 

The functions d in F are additional constraints to be 
included in the maximisation of F. In this paper, we use 
only one such constraint which is that we require the model 
luminosity to match the luminosity (L) of the galaxy being 
modelled, that is LWi — L, or more concisely Wi = 1. 
We therefore take 



Ci 



A s 



' N 



Wi 



(9) 



where A SU m is a positive parameter 

2 



Given the definitions of x > S and d and noting that 
for the purposes of this paper we do not use regularisation 
([i = 0, see section [2^3)) . G(w) from equation [2] can now be 
written 



G(w) 



K 

£ 

k 



At, — —Ak,j + A s 




(10) 



where Kk,j(ri,Vi) has been abbreviated to Kk,ji and we 
have assumed that, for all observables, 1 particle contributes 
only at 1 position j. 

Finally, model observables (and thus particle weights) 
are subject to noise as the numbers of particles contributing 
to the observables vary. This noise is suppressed by replacing 
Ak,j in G(w) by an exponentially smoothed version Ak,j 
given by 



d x , 
— A kJ = a{A ktj 



(11) 



where a is a small positive parameter. The smoothed A k j 
can be used to calculate a smoothed version y^j of the model 
observable, 



Vk,j = Yk.j +a k jAk ] j. (12) 

Assum ing the calculation suggested in ISver fc Tremainel 
is followed, the effect of exponential smoothing is as 
if the number of particles has been boosted by two or three 
orders of magnitude. 



2.3 Regularisation 

Given that the number of particles typically exceeds the 
number of observations by several orders of magnitude, de- 
termining the particle weights using the M2M method is 
an ill-posed problem with potentially many solutions. Reg- 
ularisation (S) is introduced to avoid overfitting the data 
and to reduce the number of solutions to those that satisfy 
the regularisation condition. In the M2M case, regularisation 
behaves to smooth the M2M model either at a global level 
(entropy based functions) or at a local level (second deriva- 
tive smoothing). The amount of regularisation is controlled 
by the parameter p which is used to achieve some desired 
balance between fitting the observed data and delivering a 
smooth solution. Quite what that balance should be is up to 
the individual researcher to determine. Entropy based reg- 
ularisation introduces another consideration in that it pre- 
vents the particle weights (wi) moving too far from their 
priors (m.) and thus the choice of priors has a major influ- 
ence on the model solution. 

For the investigation in this paper (determining the bar 
angle and pattern speed), we are using the M2M method as a 
'black box' and are more concerned with differences between 
models, as represented by the end of run model x 2 values, 
rather than the detailed solutions represented by individ- 
ual models. The bar angle and the pattern speed are global 
properties of the galaxy being modelled and we do not wish 
the values we determine to be influenced by the behaviour 
of any regularisation condition. We thus choose deliberately 
to model wi th no regularisation. This approach is consistent 
with that in lLong fc Maoll|2012D where we successfully used 
the M2M method as a 'black box' to determine a different 
global property (the mass-to-light ratio) of the galaxies be- 
ing modelled. 

As will be seen in the rest of this paper, the values of 
the bar angle an d patte rn speed we achieved are consistent 
with IShen et al.l ((2010) and we meet our objectives with- 
out employing regularisation. We will comment further on 
regularisation in our conclusions in section [6] 



3 MODELLING WITH MILKY WAY 
OBSERVABLES 

The two main issues to be addressed here are 

(i) modelling with Galactic field based observables where 
the fields are non-contiguous and specified in terms of Galac- 
tic longitude and latitude and can not be modelled as par- 
allel projections, and 

(ii) using a rotating frame (representing the bar) for orbit 
integration. 

Both these issues are straightforward to resolve in a M2M 
context. 



4 R. J. Long, Shude Mao, Juntai Shen and Yougang Wang 



•M««M«M$M«HMMt 



1( 25 20 15 10 5 -5 -10 -15 

I (degrees) 

Figure 1. The available BRAVA fields. There is no significance 
to the colour coding other than to distinguish adjacent or over- 
lapping fields. 



3.1 Observables 

We take the measured kinematic data for our models from 
the publicly available BRAVA dat43 (downloaded in De- 
cember 2011). The M2M observables thus comprise mean 
radial velocity and velocity dispersion squared in a number 
of Galactic fields specified by G alactic long i tude a nd lati- 
tude (I, b). For consistency with IShen et al. (|2010t ). we do 
not use the individual stellar velocity measurements within 
each field as M2M constraints. We do however include the 
b = —6° BRAVA m easurements which were not available 
to IShen et all <|20ld ). Figure [T] shows the individual fields. 
For modelling purposes, we define the boundaries of each 
field as the minimum and maximum longitude and latitude. 
In total, 86 {l,b) fields where downloaded of which 78 dis- 
tinct fields are used in our models. The difference is due 
to overlaps occuring in 15 of the fields with (I, b) values of 
(-8°, -4°), (4°, -5°), (0°,-3.5°), (0°,-6°), (10°, -6°) and 
(6°, —4°). After visual examination of the degree of overlap, 
8 of the overlapping fields were arbitrarily removed from the 
modelling. Given the low number, the fact that the spread 
of measurement points was not impacted, and the nature of 
the exercise (results consistency check only) , further investi- 
gation of the overlapping fields was considered unnecessary. 

'Measured' luminosity observ ables are calculated di- 
rectly from the IShen et al.l {2010) particle data. Luminos- 
ity density is taken to be axisymmetric (with z-axis, the 
symmetry axis) and is formed by binning particles in (R, z) 
where R is the equatorial plane projected radius. All par- 
ticles at this stage are assumed to have equal luminosity. 
The observed values are calculated out to R — 13kpc. For a 
projected luminosity observable, we use the fractional field 
luminosity where the fields match the BRAVA (I, b) fields, 
and construct it by counting particles accordingly. To take 
into account that the measured observables are determined 
by binning particle data (and are therefore subject to noise), 
and to allow the M2M models the opportunity to retain 
the bar shape (which is not axisymmetric), we set the error 
terms for the observables to a high relative error of 20%. 

It is worth noting that, at any one time, only « 2.5 x 



10 4 particles (out of ~ 10 6 particles in total) are present 
in the BRAVA fields and contributing directly to the model 
observable calculations. For binning particles into fields, we 
implement a scheme which first eliminates quickly particles 
which are not in any field before performing a more detailed 
field selection. 

The units we use within the M2M models in this paper 
are kilo-parsecs for length, 10 7 years for time, and mass in 
units of the solar mass Mq. 



3.2 Kernels 

We list the kernels used for calculating model observables. In 
the equations below, L is the total luminosity of the model, 
Vj is the volume of observable bin j, v ra diai,i is the radial 
velocity (with respect to the Sun's position) of particle i, 
and Lj is the fractional luminosity of field j. 

(i) luminosity density 

(ii) fractional field luminosity 
K Si = l 

(iii) 1st radial velocity moment 

^radial, i 



Kn = 



Li 



(iv) 2nd radial velocity moment 



K u = 



(13) 



(14) 



(15) 



(16) 



We calculate the model radial velocity dispersion 
squared directly from the 1st and 2nd velocity moments us- 
ing 



« 2 j - (vj) 



(17) 



where, for bin j, <jj is the dispersion, v 2 j is the 2nd moment, 



3.3 Gravitational Potential 

As in IShen et all (|201Ch . the gravitational potential com- 
prises 2 components, a dark matter halo given by 



r 

1 + m 



(18) 



http://brava.astro.ucla.edu/data.htm 



where V c = 248 km/s and R c — 15.2 kpc, and a potential 
derived from the luminous matter particles. 

We do not u se self-gravitation of the particles as in 
IShen et all \20ldh and comment on this later in our conclu- 
sions (section [BJ. Instead, we calculate the luminous matter 
potential and its associated accelerations once from the par- 
ticle positions using direct summation and hold the values 
in interpolation tables. The tables use cylindrical polar co- 
ordinates with a pseudo-logarithmic radial coordinate. We 
use tri-linear interpolation to determine the function val- 
ues as required at specific po sitions. We initially used the 
self-consistent field method of iHernquist fc Ostrikerl (l992) 
to determine the potential from the particles but stopped 
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when we discovered that, under certain circumstances, den- 
sities calculated via the method were negative. Why this 
occurred is currently unresolved. 

Direct summation is straightforward to parallelise using 
a graphics processor unit (GPU) and we currently achieve 
a four-fold increase in performance over a single computer 
processor. With experience, we would expect this increase 
to become greater. 

The total mass of the particle model is 4.25 x 1O 1O M0. 
For implementation reasons (the particle model is a mass 
model and the M2M theory in section [2] is framed in terms 
of luminosity), the M2M models assume a luminous matter 
mass-to-light ratio of 1. 



3.4 Particle initial conditions 

The initial positions and veloc ities for the partic les are the 
end of run values taken from IShen et al l (|20ld ). Plots of 
the projected initial positions of the particles are shown in 
Figure [2] In these plots, the bar lies along the x-axis. During 
M2M modelling, the particles are first rotated to the desired 
bar angle before orbit integration starts. 

The initial weights of the particles are set as 1/N where 
N is the number of particles. 



3.5 Rotating frame kinematics 

It is convenient to think of the M2M model in terms of 3 
coordinate frames (inertial, rotating and constraints) with a 
common origin at the Galactic centre. 

The inertial frame represents the overall model with 
the Sun stationary on the positive x-axis at 8.5 kpc from 
the Galactic centre. The x-y plane is the equivalent of the 
Galactic equat orial plane a nd is shown diagrammatically in 
Figure [3] The IShen et all |2010T l particle data has the bar 
positioned along the x-axis initially with a clockwise sense of 
rotation. For modelling purposes, the bar is first positioned 
at some desired bar angle prior to orbit integration starting. 

The rotating frame rotates at a given pattern speed 
(f2 p ), fixed per modelling run, and contains the bar (and 
indeed all the particles). Orbit integration takes place in 
the rotating frame and u ses a leapfrog scheme based on 
iPfenniger fc Friedlil (|l993l ). We use the leapfrog scheme 
in 'drift-kick-drift' form with a constant time step of 
2 x 10~ 3 time uni ts . Jac obi's integral is conserved (see 
iBinnev fc Tremaind l|200&t )) to a relative accuracy of ~ 
10 -2 . This value is the mean maximum relative accuracy 
taken across all particles for the duration of a modelling run. 
Given we are using interpolation tables to hold the potential 
and its accel erations, and are o nly seeking consistency with 
the results in IShen et al. (2010), the value is satisfactory for 
our purposes. We orbit the particles for 125 time units. For 
Q p = —40 km/s/kpc, this equates to ~ 8 full revolutions of 
the bar. 

Observable constraints must only be applied with the 
bar at the desired bar angle to the Sun to Galactic centre 
line. Conceptually, the constraints application frame must 
therefore be rotated as the bar rotates to keep the bar angle 
at the desired value. This is the equivalent of either rotating 
the Sun's position or reversing the bar's rotation for the 
purposes of applying the constraints. In practice, to avoid 



rotation 




Figure 3. Schematic showing the inertial equatorial plane of the 
M2M model. $bar is the bar angle and has a negative value when 
mimicking the Galactic bar. Given the sense of rotation of the 
bar is clockwise, the pattern speed (f2 p ) is also negative. Positive 
Galactic longitudes are on the left of the schematic. 



Table 1. M2M parameters 



Parameter 


Value 


Related constraint 


c 


2.5 x 10 -3 




a 


5.0 x 10~ 2 











^sum 


10 4 


sum of particle weights 




10" 3 


luminosity density 


Afl 


10" 3 


fractional field luminosity 


Av 


10~ 4 


radial velocity 


Avd 


10" 4 


radial velocity dispersion 



The M2M parameters and their values. The same values are 
used in all modelling ru ns. The rationale for the observable A 
•arameters is recorded in lLong fc Maol ll2010h and lLone fc Maol 
20121 '). A sum is set by experimentation to ensure 

J2m = i-oo. 



manipulating the constraints or their fields, we choose the 
latter. 



3.6 Parameter setting 

In Table [1] we list all the parameters identified in section 
12.21 and the values we use in our M2M models. The values 
are the same for all models. The rati onale for the Xk pa - 
ram eters in eq uation [3] is recorded in lLong fc Maol (|2010) 
and lLong fc Mad (|2012T ) and is not repeated here. These 
parameters have only been set to an order of magnitude. 
No fine tuning, or sensitivity analysis, has been undertaken. 
A sum has been determined by experimentation to ensure 
^2 Wi = 1.00 for all models. 



4 RESULTS 

By varying the bar angle #bar between —10° and —70° with a 
10° interval, and the pattern speed Q p between km/s/kpc 
and —70 km/s/kpc with a 10 km/s/kpc interval, we con- 
struct a set of 56 M2M models. Using these models, we cre- 
ate a model x 2 surface where x 2 i s as m equation [3] and 

and 
-30° 



involves the Xk values. By marginalising \ over $ba 
rip individually, we obtain the best fit values of i9b ar ~ 
and Q v 



-40 km/s/kpc. These values are consistent with 
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Figure 2. Logarithmic fractional projected particle counts showing the initial spatial distribution of the particles. The bar lies along the 
x-axis (indicated by the red line in the left hand panel). During modelling, the particles are first rotated to the desired bar angle before 
orbit integration starts. 



the values IShen et al l (|2010h determined for their N-body 
model. If we now take only the kinematic components from 
X 2 , without their corresponding At, values, and repeat the 
marginalisation exercise, we arrive at the same values for 
#bar and ftp. The \ 2 surfaces and marginalised plots are 
shown in Figure 3] It is clear from the figure that the bar 
angle is only weakly determined by the kinematic observa- 
tions alone, and is better determined when the luminosity 
constraints are included. 



The best fit model (0 b 



-30° and fi p 



— 40km/s/kpc) reproductions of the BRAVA field mean ra- 
dial velocities and veloci ty dispersions are given in Figure [5] 
The equivalent figure in IShen et al l l|201Ch is Figure 2. The 
model reproductions of the kinematics are quite smooth and 
do not suggest that regularisation should have been used. 
The \ 2 P er degrees of freedom values for the best fit model 
are 2.55 for the radial velocities and 1.09 for the velocity dis- 
persions. It may be seen from Figure [S] (see the b = —6° and 
b = —8° panels) that the scatter in the field radial velocity 
values and the size of the error bars is causing the M2M 
models to be unable to reproduce all the measured veloc- 
ity values, resulting in higher \ 2 values not just in the best 
fit model but in all models. The \ 2 per degrees of freedom 
values for luminosity density and fractional field luminosity 
are both low at « 0.1. This is almost certainly due to using 
high relative errors (20%) as noted in section ETT1 However, 
given the consistency objective in section [T] had been met, 
no experimentation took place with lower error values. For 
comparison with Figure Figure [5] shows the variation in 
observable reproduction for pattern speeds either side of the 
best fit value, with the bar angle at its best fit value. 

Notwithstanding that an axisymmetric distribution has 
been used in constructing the luminosity density constraint, 
the bar structure has been maintained during orbit integra- 
tion as can be seen from Figure [7] Weight convergence for 
our best fit model is high with 99% of particles having con- 
ve rged weights (the d efinition of weight convergence is given 
in lLong fc MacfeoioT ). The final weights of the best fit model 
particles remain close to their initial values (1.02 x 10 -6 ) and 
are in the range 5.9 x 10 -7 < Wi < 2.8 x 10~ 6 . The variation 
in weight range is most noticeable along the pattern speed 
axis where the minimum of the range is at a maximum for 
the best fit pattern speed Q p w — 40km/s/kpc and decreases 
on either side to w 10 -7 at — 70km/s/kpc and « 10 -8 at 
km/s/kpc. As a final comment, the narrow weight range 




n: 



-10 10 

X (kpc) 



Figure 7. Projection onto the equatorial plane of the logarithmic 
particle weight distribution for our best fit model ($bar ~ —30° 
and Sip pa —40 km/s/kpc). The bar structure (see Figure [2]l is 
maintained during orbit integration and is in the expected posi- 
tion at end of run (indicated by the red line). 



of the best fit model is not unexpected given the way the 
luminosity data was constructed and the use of the BRAVA 
data in the N-body model. The variation with pattern speed 
we take as a positive indicator that the M2M method used 
as a 'black box' is able to discriminate successfully between 
top level parameter values. 

As a guide to the computing resources required, one 
M2M modelling run, using 32 processing cores, takes ~ 60 
minutes to complete. Assuming the M2M models in this ex- 
ercise are run sequentially, the full set takes 56 elapsed hours 
to complete. We did not encounter t he da ta packet fragmen- 
tation issue noted in lLong fc Mar] l|2012l ) as the number of 
measurement points (equal to the number of fields) is low. 
Whether other M2M modellers will hit the fragmentation is- 
sue depends on their M2M implementation and its use of the 
Message Passing Interface (MPI), the variant of MPI used, 
and the network configuration of the hosting computer sys- 
tem, as well as the number of observable measurement points 
for the stellar system being modelled. 
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Figure 4. x 2 surface plots and marginalised plots for the M2M models using both density and kinematic constraints (left hand column), 
and for the kinematic constraints only (right hand column). In the left hand column, x 2 is calculated using equation [3] involving the 
parameters. No X k values are used in calculating the kinematic x 2 m the right hand column. In the marginalised plots, the data points 
are shown in red; the blue lines are for visual guidance only. tt p is the pattern speed and 6>bari the bar angle. 



5 DISCUSSION ON COMBINED N-BODY 
AND M2M MODELLING 



We now consider further combining N-body and M2M mod- 
elling. Our initial thoughts on starting the current exercise 
were that we could forsee a combined modelling process with 
repeating interleaved N-body and M2M phases with each 
phase seeking to refine the preceding phase. This paper rep- 
resents the first step in developing that process in that it 
seeks to ensure that M2M modelling delivers results consis- 
tent with an established N-body stellar dynamical model. 



However the interleaving approach is not the only way for- 
ward and we now describe some alternatives. 
There appear to be four main options: 

(i) keep N-body and M2M modelling completely separate 
so that they can be used to perform an independent check 
on each other, 

(ii) loosely couple the 2 methods so that the output from 
an N-body model can be fed into a M2M model (with the 
ability to iterate between model types), 

(iii) absorb N-body modelling into M2M modelling, and 

(iv) absorb M2M modelling into N-body modelling. 
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Figure 5. Reproduction by the best fit M2M model (6>bar ~ —30° and Q p fa —40 km/s/kpc) of the observed velocities and velocity 
dispersions along the minor axis (1 = 0°) and by latit ude (ft = —4°, —6°, —8°). The black solid lines represent the model and the observed 
values arc in blue. Sec Figure 2 in IShen et al.l feoiCft for a comparison with the N-body results. Note that positive Galactic longitude is 
to the right in the plots to aid that comparison. 



These options are not necessarily disjoint. Practically, in a 
software engineering sense, the last three options may not 
actually turn out to be that dissimilar. For a given mod- 
elling scenario it may be appropriate to reuse one model 
type with multiple models of the second type, as in the cur- 
rent paper where one N-body model was used with many 
M2M models. Clearly the N-body model could have been 
re-established for every M2M model but that would not of- 
fer a computationally cost effective solution. In other words, 
even in a combined model, there is a case for keeping the 
N-body and M2M modelling phases separate. 

The one-to-many scenario is appropriate when the rea- 
son for modelling is the determination of some global prop- 
erty or attribute of a stellar system, for example, the mass- 
to-light ratio or the bar pattern speed. In this case, individ- 
ual models are less important than the comparison between 
the models. However, there are cases where the individual 
models take precedence (for example, in investigating or- 
bital structures or velocity dispersion anisotropy), and the 
relationship between the phases comes closer to one-to-one. 
The case then for keeping the phases separate appears to 
reduce. It must be remembered that the M2M method can 
only weight existing, pre-specified orbits / particles. It does 
not as yet include a mechanism to modify the orbits dy- 
namically. Currently the only way to do so is to change the 
gravitational potential or the initial conditions of the parti- 
cles by starting a new model. In the context of the current 
discussion, this means invoking the N-body phase. 

As with all categorisations, situations arise which do not 
fit completely into any one category. Within a M2M model, 



orbiting the particles under self-gravitation perhaps with an 
additional (dark matter) potential might be appropriate. For 
example, in the currrent investigation, if we had had the ca- 
pability within our M2M implementation, it would have re- 
moved the need to construct interpolation tables. Whether 
the self-gravitation of particles is regarded as an integral 
part of M2M modelling or the absorption of N-body mod- 
elling into M2M modelling is open to debate. 

For the fourth option, conceptually at least, without the 
benefit of practical experience, adding a M2M capability to 
an existing N-body software package would appear to be 
straightforward. Extrapolating, it should also be straight- 
forward to add a capability to an existing N-body / gas 
dynamics code, that is to include a M2M style of handing 
observables. 

Regardless of the approach adopted, there is still the 
need to determine how the two types of modelling should 
be interfaced. Certainly initially, this is seen as a manual 
process requiring assessment of the modelling so far and a 
judgement as to whether or not the next phase should be 
invoked. Automation becomes possible once the 'rules of en- 
gagement' become clear. In addition, quite how time is to 
be regarded must be addressed. In N-body models, it is not 
unusual for the evolution of the model with (astrophysical) 
time to be important, for example in understanding how 
galactic bars are created and evolve. In a M2M model, time 
is no more than part of the mechanism used in adapting the 
particle weights so that observations of a real galaxy, taken 
at some specified astrophysical time, are reproduced. Also, 
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Figure 6. Reproduction by the best fit M2M model (6>bar ~ —30° and Q p fa —40 km/s/kpc) of the observed velocities and velocity 
dispersions along the minor axis (/ = 0°) and by latitude (ft = —4°, —6°, —8°). The black solid lines represent the best fit model and the 
black points with error bars, the observed values. With the bar angle at the best fit value, the blue and red lines show the observable 
reproduction for patterns speeds of Q p = km/s/kpc and Q p = —70 km/s/kpc respectively. 



dynamical stability over time of the resulting models will 
need to be considered. 

To conclude this section, we have identified some alter- 
ative approaches for combining N-body and M2M modelling, 
and pointed out some of the issues which will need to be ad- 
dressed in future investigations. 



6 CONCLUSIONS 

We have met the objectives we set out in the introduction, 
section [1] We have extended the made-to-measure method 
to handle Galactic field observations involving non-parallel 
projections, and used it with rotating frame kinematics to 



model the Galactic bar. Utilising the lShen et al.l (|2010h N- 
body particle model of the Galactic bar and bulge, and kine- 
matic data from the BRAVA survey, we have used the made- 
to-measure method to determine the best fit bar angle and 
patt ern speed m a tching the data. Our results are consistent 
with IShen et all (|2010r i. Whilst this is pleasing, it was not 
wholly unexpected given the construction of the M2M mod- 
els. The rotating frame orbit integration is not specific to 
the Milky Way and co uld be used with any rotating galaxy. 
The lShen et al.l l|2010h bar pattern speed of w 40 km/s/kpc, 
corroborated by the M2M models in this paper, is lower 
than the usually quoted value of 60 km/s/kpc for the Milky 
Way. It is unclear whether such a lower pattern speed is 
consistent with differ ent observations such as gas dynamics 
( Biss antz et all l2003h and local kinematics (|Dehnenl l2000; 
iLiu et all 120121 ). The IShen et all (|2010r i pattern speed is 



however close to the valu e (42 km/s/kpc) determined by 
IWeiner fc Sellwoodl (|l999h using gas kinematics, as is the 
bar angle (34°). 

Regularisation in the M2M method was not used and 
from our modelling results was not required. This does not 
mean that we do not expect to use regularisation in future 
modelling. We will decide based on the requirements of the 
investigation to be performed and the results we are achiev- 
ing. Our implementation of the M2M method would benefit 
from having a particle self-gravitation capability to avoid 
the need to take a 'snap shot' potential from an N-body 
system, and to avoid the use of interpolation tables. We will 
address this in the near future. 

Future simple extensions to the M2M method in- 
clude modelling the individual stellar velocity measure- 
ments, rather than field averages, by enhancing the theory in 
lLong fc Maol (|201Ct ) to handle a dark matter contribution to 
the line-of-sight velocity distribution. It would be interest- 
ing to understand the impact, if any, on the results of mod- 
elling with field averages and individual measurements both 
in combination and separately. A related extension would be 
to include proper motion data, for example fro m the Opti- 
cal G ravitational Lensing Experiment (OGLE) (|Sumi et al.l 
l2004h . in M2M models of the Galaxy. Such extensions will 
be valuable once data from the impending GAIA satellite 
become available. It is not necessary however to wait for 
GAIA before utilising the extensions to produce improved 
models of the Milky Way. 

We have presented a brief analysis of the options and is- 
sues involved in combining N-body and M2M modelling, and 
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have taken a first step in that direction. Of the issues, the 
dynamical stability of the resulting model may prove to be 
the m ost challenging. For example, the model of lWang et al.l 
|2012t ) is able to reproduce the BRAVA kinematic observa- 
tions but appears to be unstable on timescales longer than 
1 Gyr. 



ACKNOWLEDGEMENTS 

The authors wish to thank the anonymous reviewer for the 
prompt production of their very helpful comments. Com- 
puter runs were performed on the Laohu high performance 
computer cluster of the National Astronomical Observato- 
ries, Chinese Academy of Sciences (NAOC). RJL and SM ac- 
knowledge the financial support of the Chinese Academy of 
Sciences and NAOC. The research presented here is also par- 
tially supported by the National Natural Science Foundation 
of China under grant numbers 11073037 (JS), Y011061001 
and Y122071001 (YGW) and by the 973 Program of China 
under grant number 2009CB824800 (JS). 



Rich R. M., Reitzel D. B., Howard C. D., Zhao H., 2007, 

ApJL, 658, L29 
Sellwood J. A., 2012, ApJ, 751, 44 

Sellwood J. A., Wilkinson A., 1993, Reports on Progress 

in Physics, 56, 173 
Shen J., Rich R. M., Kormendy J., Howard C. D., De Pro- 

pris R., Kunder A., 2010, ApJL, 720, L72 
Sigurdsson S., Hernquist L., Quinlan G. D., 1995, ApJ, 446, 

75 

Sumi T., Wu X., Udalski A., Szymahski M., Kubiak M., 
Pietrzyriski G., Soszyhski I., Wozniak P., Zebruh K., 
Szewczyk O., Wyrzykowski L., 2004, MNRAS, 348, 1439 
Syer D., Tremaine S., 1996, MNRAS, 282, 223 
Wang Y., Zhao H., Mao S., Rich R. M., 2012, ArXiv e- 
prints 

Weiner B. J., Sellwood J. A., 1999, ApJ, 524, 112 



REFERENCES 

Binney J., Tremaine S., 2008, Galactic Dynamics: Second 

Edition. Princeton University Press 
Bissantz N., Debattista V. P., Gerhard O., 2004, ApJL, 

601, L155 

Bissantz N., Englmaier P., Gerhard O., 2003, MNRAS, 340, 
949 

Das P., Gerhard O., Mendez R. H., Teodorescu A. M., de 
Lorenzi F., 2011, MNRAS, 415, 1244 

de Lorenzi F., Debattista V. P., Gerhard O., Sambhus N., 
2007, MNRAS, 376, 71 

de Lorenzi F., Gerhard O., Saglia R. P., Sambhus N., De- 
battista V. P., Pannella M., Mendez R. H., 2008, MNRAS, 
385, 1729 

Dehnen W., 2000, AJ, 119, 800 

Dehnen W., 2009, MNRAS, 395, 1079 

Fux R., 1999, A&A, 345, 787 

Hernquist L., 1993, ApJS, 86, 389 

Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375 

Howard C. D., Rich R. M., Reitzel D. B., Koch A., De 
Propris R., Zhao H., 2008, ApJ, 688, 1060 

Jourdeuil E., Emsellem E., 2007, in Kissler-Patig M., Walsh 
J. R., Roth M. M., eds, Science Perspectives for 3D Spec- 
troscopy Scalable N-body Code for the Modeling of Early- 
type Galaxies, pp 99-103 

Kunder A., Koch A., Rich R. M., de Propris R., Howard 
C. D., Stubbs S. A., Johnson C. L, Shen J., Wang Y., 
Robin A. C, Kormendy J., Soto M., Frinchaboy P., Re- 
itzel D. B., Zhao H., Origlia L., 2012, AJ, 143, 57 

Liu C, Xue X., Fang M., van de Ven G., Wu Y., Smith 
M. C, Carrell K., 2012, ApJL, 753, L24 

Long R. J., Mao S., 2010, MNRAS, 405, 301 

Long R. J., Mao S., 2012, MNRAS, 421, 2580 

Morganti L., Gerhard O., 2012, MNRAS, 422, 1571 

Pfenniger D., Friedli D., 1993, A&A, 270, 561 

Rattenbury N. J., Mao S., Debattista V. P., Sumi T., Ger- 
hard O., de Lorenzi F., 2007, MNRAS, 378, 1165 



